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Abstract 

We derive an analytical expression for the effective force between a pair of macrospheres 



i 

immersed in a sea of microspheres, in the case where the interaction between the two unlike 

o 

species is assumed to be a square well or a square shoulder of given range and depth (or height). 



This formula extends a similar one developed in the case of hard core interactions only. Qual- 
itative features of such effective force and the resulting phase diagram are then analyzed in the 
limit of no interaction between the small particles. Approximate force profiles are then ob- 



tained by means of integral equation theories (PY and HNC) combined with the superposition 



approximation and compared with exact ones from direct Monte Carlo simulations. 



^ 1 Introduction 

Colloidal systems are ubiquitous and have attracted a growing interest in the last decades. The 
impact on everyday life is vast and ranges from food 1 to materials, from biology to photonic 
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crystals. 2 On a more fundamental level, experiments on colloids have been the testing ground for 
several physical phenomena, from thermodynamics 3 ' 4 to the glass transition. 5 In general, the in- 
terest comes from the possibility of tuning macroscopic properties by encoding them at the level 
of the interactions between the constituents of colloidal solutions. For these reasons these systems 
have dominated the scientific discussions in these fields during the last years. In particular, the 
possibility of engineering systems with tunable microscopic properties allowed the observation of 
exotic phenomena like the Alder transition for hard spheres, 6 the metastable liquid-liquid phase 
separation, 7 the existence of two distinct glassy phases 8 ' 9 for short ranged attractive potentials and 
the emergence of a thermodynamically stable cluster phase. 10 ' U 

As most soft-matter systems, colloids are characterized by a large number of degrees of freedom 
spanning several length and time scales. In several circumstances, however, it is possible to inte- 
grate out some of these degrees of freedom to gain more insight. In fact, colloids are one of the 
prototypical systems where a coarse-grained approach can be extremely fruitful. Typically, the 
extra degrees of freedom are mapped into some effective interaction among colloidal particles, 12 
enabling the use of the full arsenal of statistical mechanics. In particular, most of the results for 
simple liquids (integral equations, perturbation theory, mode coupling theory, etc.) can be success- 
fully used to address questions of great relevance for colloidal systems. 

A classical example of such success is the case of depletion interactions. In the case of a binary 
mixture of colloidal hard spheres, when the diameters of the two species are very different, an 
entropic force starts to set in. 13 This results in a net attraction between the colloids belonging to 
the largest species. Asakura and Oosawa 14 and, independently, Vrij 15 derived an effective interac- 
tion potential casting the problem of a binary mixture into an effective single-component system. 
This very simple result was the beginning of a series of investigations (see for instance (16,17)). 
It was possible, for the first time, to investigate a system of particles interacting with an attraction 
whose range could be tuned, something that is not feasible for atomic or molecular systems. The 
profound effects both on thermodynamics and the dynamics of a short range attraction are now 
well established. 



The use of effective interactions, far from being restricted to colloidal systems, has important im- 
plications also for proteins, especially for what concerns crystal nucleation and phase behavior. By 
modeling proteins as short range attractive colloids it has been possible, for example, to rationalize 
the enhancement of crystal nucleation in the proximity of the liquid-liquid critical point. 18 This 
extremely simplified modeling has been also useful in the case of binary mixtures of eye-lens pro- 
teins, whose experimental behavior was successfully modeled using a simple hard sphere potential 
with specific square well attractions. 19 ' 20 Here, a fundamental role was played by the intensity of 
the interspecies attraction: by varying it, it was possible to reduce the instability with respect to 
demixing due to depletion. Apart from the medical and biological implications, this work showed 
how a modification of the mutual interaction among the components can alter the phenomenology 
of the whole mixture. 

The aim of this paper is to rationalize this fact in terms of effective interactions. We investigated 
the effective interactions between the larger component of an asymmetric binary mixture of hard 
spheres where the potential between the two components has the form of a square well or a square 
shoulder. To derive our results we followed a reasoning similar to Attard's derivation of the effec- 
tive force between two hard spheres immersed in a sea of smaller ones. 21 We generalized Attard's 
argument to the case of square well (or shoulder) interspecies interaction. The calculation of the 
actual force, as in (21), requires the evaluation of the density profile of the smaller spheres around 
a pair of larger ones. We performed this step following different methods. We shall first propose 
a simple approximation that treats the smaller component as an ideal gas. In this way, we shall 
find a generalization of the Asakura-Oosawa (AO) analytical expression of the force to the case in 
which the mutual interaction potential is an attractive square well or a repulsive shoulder. Within 
such approximation we'll explore, by means of first order thermodynamic perturbation theory, the 
qualitative changes in the phase diagram determined by the mutual interaction. The analogy with 
the case of binary mixtures of eye-lens proteins are evident. As for the purely repulsive case, the 
approximation above holds only when the interactions between the smaller particles are negligi- 
ble. Thus we improved the results of the force calculation using an integral equation approach 



combined with the superposition approximation for the case of hard sphere interaction between 
the smaller species. All the aforementioned theoretical results have been tested with respect to the 
"exact values" of the force as calculated directly from new Monte Carlo simulations. 
The paper is organized as follows: in Section 2 we present our generalization of Attard's deriva- 
tion of the effective force in the case of square well/shoulder interactions between the species of 
a binary mixture; in Section 3 the methods used to estimate the density profile needed to com- 
pute the effective force are described. Details about the parameters used are also provided and a 
description of the qualitative effect of the effective interactions on the phase diagram is given; in 
Section 4 quantitative results obtained with the different methods are shown and discussed. Finally 
in Section 5 conclusions of this work are drawn. 

2 Force in the square well/shoulder case 

Consider an asymmetric binary mixture (we label with 1 the larger component and with 2 the other 
one) and let the interaction potential between two particles of the larger species 0n(r) and the 
interspecies potential <j>i2(r) have a square-well or square shoulder form, that is 



<M r ) = < 



°o if r < Oij 

£ij if Oij < r < Xij Oij , (1) 

if Xij Oij < r 



where Oij denotes the range of the hard core interaction, £ ;J is a positive (negative) energy mea- 
suring the depth (height) of the well (shoulder) and Ay determines the width of the well (shoulder) 
with respect to the range of the hard core. We fix £n = and leave the interaction potential of 
the smaller species 022 (r) unspecified. We now wish to calculate the effective force between two 
particles of type 1 immersed in a sea of particles of type 2. 

The problem in the particular case of £12 = (the hard sphere limit) has already been studied in 
detail 21 ' 22 in the past. If two spheres of type 1 are fixed in (the origin) and R respectively, the 



value / of the radial component of the effective force between each other can be written as (21): 

f{R) = -—o 2 n p(S ai2 ;R)cos0sm0d0, (2) 

p Jo 

where /3 is the inverse temperature l/kgT, S ai2 is a vector of length 0\2 whose tail is in and 
6 is such that S ai2 R = R 0\2 cos 9 (the meaning of the parameters that appear in Eq. (2) is further 
elucidated in Appendix A and Figure 1). Thus the force depends on the density p(r) of small 
particles around the pair of macrospheres. 

In Appendix A a derivation similar to Attard's is used to obtain the expression of the effective 
force between two macrospheres when the interaction between the species 1 and 2 is of the kind 
described by Eq. (1): 

f(R) = — g- lo 2 2 J p(S ai2 ;R)cos0sm0d0 + 

(l-e^)X 2 af 2 J%(S Xan ;R)cosesmdde) , (3) 

where the shorthand notation A a\% = ^n^n has been used, S^ au is defined in a fashion similar 
to S 0n and a representation of all the parameters is given in Figure 1. 

Our expression differs from the hard sphere case ( Eq. (2)) because of the presence of a second 
term. This is due to the introduction of the square well (shoulder) of finite depth (height) \S\2\ and 
has the same form of the previous term but a weight 1 — e^ £n , which is negative for square well 
potentials. This dependence of the force on the temperature is a by-product of the non-zero value 
of £12, which introduces a natural energy scale that is absent in the athermal hard sphere system. 
Thus, the force is not entirely entropic as in Eq. (2) but has also an energetic origin. The expression 
of the force in Eq. (3) correctly reduces (as it should) to Attard's formula in the £12 = (hard 
sphere) limit and to the same expression and a larger a = X<5\2 in the limit £12 — > — °° (infinitely 
high shoulder). We emphasize the fact that Eq. (3) yields the exact force if the correct form of 
the density is known regardless the interaction potential between the small particles. In this paper 



we shall restrict ourselves to the two cases in which the small particles are either non-interacting 
(022 (r) = 0) or hard spheres (so that (foii?) will take the form of Eq. (1) with £22 = 0). 

3 Determination of the density and force 

In order to determine the effective force the density p(r) was obtained using three different meth- 
ods: the dilute gas approximation, integral equations together with the superposition approximation 
and Metropolis Monte Carlo simulations. 

3.1 Dilute limit (DL) approximation 

If species 2 is dilute we can approximate the density around a pair of macrospheres with that of an 
ideal gas: 

p Dh (r;R)=pe- pv{r;R \ (4) 

where V(r;R) is the potential energy of a particle of type 2 centered in r due to the presence of the 
fixed pair of type 1 (see also Eq. (21)). 

This is equivalent to the Asakura-Oosawa approximation in the case £12 = 0. We expect such 
approximation to be less accurate in more coupled regimes but still able to give qualitative infor- 
mation about the force profile. In App. Appendix B we make use of Eq. (4) and obtain an analytic 
expression of the effective force. For a wide choice of the parameters £12, C12, kon it has features 
like those appearing in Figure 2. The profile obtained when £12 = is exactly the same short range 
depletion attraction given by the Asakura-Oosawa approximation. 21 The effect of negative values 
of £12 is that of increasing both the range and the strength of the effective force, approaching the 
g 12 = —00 solution mentioned above. In this case the qualitative features are essentially the same as 
in the £12 = case. Positive values of £12 determine a completely different behavior: increasingly 
large values determine the onset of a repulsion at short distances (maximum repulsion occurring 
at R = 2o\2) and of a strong attractive force at intermediate distances (with maximum attraction 
found at R = <3\2 -\-Xo\2). 



The physical origin of such form of the effective interaction can be understood (at least quali- 
tatively) by means of the argument below. 

For sake of simplicity, let's consider the £12 = case first. In this case one macroparticle and a mi- 
croparticle repel each other when they come into contact, i.e. when they are at a distance r = G\i. It 
follows that if a macrosphere is inserted in a sea of small ones, it will experience on average a force 
depending on the density of small spheres at distance r = <j\2. If a macrosphere is isolated from 
others the density distribution of small ones will be spherically symmetric so it won't experience 
any net force. If two macrospheres come close together, however, the hemispheres at r = 012 fac- 
ing each other become depleted of microspheres: the "push" from the two external hemispheres is 
not balanced anymore and the result is an effective attraction. The situation is depicted in Figure 3. 

Now consider the case when a square well or shoulder is also present. We have to take into 
account that a macroparticle and a microparticle experience an attraction or a repulsion depending 
on the sign of£yi when their separation is A <J\2. Thus, the effect of a well is to make a macrosphere 
being attracted by the density of small spheres localized at a distance r = Xa\2, while a shoulder 
will make the macrosphere be pushed away from it. Taking this into account we can understand 
the effect of £12 on /dl Ce- 
lt's easy to realize (see also Figure 3) that the case £12 < (shoulder) is similar to that where 
£12 = 0, with the novel contribution due to the density localized at the outer rim of the shoulder 
and a weaker contribution from the hard wall due to lower density of small particles inside the 
shoulder. 

The situation is more complicated if £12 > (well) and is schematized in Figure 4. If <j\2 + 
Xo\2 < R < 2Xo\2 part of the outer rim of the well is inside the well of the neighboring big 
particle and there is a strong "pull" due to the high density of small particles here (green arrows in 
Figure 4(b)). The result is an effective attraction. At lower R such attraction is counterbalanced by 
a strong "push" due to the fact that the density at the hard wall is higher where the two cores face 
each other (see blue arrows in Figure 4(c)). At the same time the attractive "pull" at the outer rim 
of the well (green arrows on the right in Figure 4(c)) is weaker because part of it is in the region not 



accessible to the small spheres (p = 0). For these reasons at small R and big enough £12 repulsion 
dominates over attraction ( Figure 4(d)). 

Qualitative dependence of the phase diagram from £12 

We studied the effect of the presence of the smaller species on the phase diagram of the larger 
component, to see if the stability of the system can be tuned by means of varying the parameter £12. 
Effects of the mutual attraction on the stability of this class of systems has previously been reported 
in the literature. 19 ' 20 In the semi-grand canonical ensemble we can write for the thermodynamic 
potential F(N,V,Z2) 

where N is the number of large particles, V is the volume of the system, zi is the fugacity of the 
smaller species, 23 Tri denotes the integration operator over the coordinates of the particles of the 



larger species / dR\ . . .dR^, A\ = h/ ' yj1itm\l '/3 its thermal wavelength and // e ff is the effective 
potential. // e ff in turn can be written as 

H eS = Hn + &, (6) 

where H\\ is the sum of the pair interactions between the larger particles T^Lj <j>n(R) an d Q. is the 
grand potential of the smaller species in a fixed configuration of the large particles. It can be shown 
that Q. can be expressed as a sum of n-body terms 12 ' 23 

a = £^„ (7) 

n 

and it is straightforward (see Appendix C) to prove that the 0- and 1-body terms are linear in the 
density of the effective component p\ and therefore do not alter the phase behavior 12 ' 23 in the case 
examined here. In what follows the 2-body term is assumed to be given by 



N 

^2 = £V DL ( J R !J ), (8) 



where 



V DL (R) = -j fMR)dR + c, (9) 

where c is chosen to guarantee that linv^oo VdlM = holds. 

First order thermodynamic perturbation theory 24 can thus be employed to explore the phase 
behavior of the effective component, in a way similar to that applied by Gast and coworkers in the 
AO case (£12 = 0). 4 We thus write the free energy per particle of the effective component F eff /N 
with 

PF eff /N = PF H s/N+^Jv DL (r)gm(r)4nr 2 dr, (10) 

where Fus/N and gns are respectively the free energy per particle and the pair distribution 
function associated to the hard sphere reference system. The former is approximated here by 
means of the Carnahan-Starling 24 expression and the Verlet-Weis 25 description is used for the 
latter. A check of the validity of the first-order approach can be performed by verifying that the 
Barker-Henderson second-order correction to F^/N 4,26 is reasonably smaller than the second 
term on the RHS of Eq. (10). The chemical potential ji and the pressure p can be expressed with 



/3 M = ^(Pi/3F eff /AO, (ID 

dpi 



Pp = Pl Pll-piPF eff /N. (12) 

Coexistence lines are found by plotting the value of the volume fraction r\\ of the bigger spheres 
at which the parametric curve in the (p,/l) plane self-intersects for various values of the volume 
fraction r\ of the smaller particles. 



3.2 Integral equations and superposition approximation (PY + S and HNC 
+ S) methods 

Another approximate way to determine the density is using integral equation theory for a fluid 
mixture in order to get the unlike (big-small) pair distribution gn(r). This in turn can be plugged 
into the superposition approximation for the density profile around the pair of macrospheres 

Ps(r;R) = gi 2 (r)gi 2 (\r-R\)p. (13) 

Note how the superposition approximation amounts to saying that the density around the pair is 
equal to the product of the densities that the particles would have around themselves if they were 
isolated. 

The pair distribution function gi2 in Eq. (13) can be determined solving the two-component 
Ornstein-Zernike equation 

2 , 

hij{r) = Cij (r) + £ / p k c ik {\r-r>\)h kj {r>)dr>, (14) 

k=\ J 

and the closure equation 

g t j(r) = e -P^j( r ) +h u( r )- c ij( r )- E ij( r ) ? (15) 

where i, j, k G {1,2}, hij,Cij,E{j(r) are respectively the indirect correlations, direct correlations 
and bridge functions, 27 p k is the density of the species k and Y*kPk = P- I n order to solve the 
integral equations we implemented a version of Gillan algorithm. 28 Special care was taken to treat 
discontinuous potentials of the form of Eq. (1). 
The well-known Percus-Yevick closure (PY) 

Eij(r) = ln[l + hij(r) - Cij {r)] - h^r) + Cij (r) (16) 
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and the hypernetted-chain closure (HNC) 

Eij(r)=0 (17) 

were chosen. 

3.3 Monte Carlo (MC) method 

The exact density can be sampled using Metropolis Monte Carlo simulations. The force can then be 
obtained with the method used in the two-sphere studies in (22) (a possibly more efficient alterna- 
tive method is that described in (29)). We briefly summarize such method here. Two macrospheres 
of diameter an are inserted at fixed positions (0,0,0) and (R, 0,0) in a cell whose dimensions are 
H along the x direction and L in the y and z directions. The same cell contains N smaller particles 
which interact with each other with a potential of the form of Eq. (1) such that £22 = 0. H and L are 
chosen so that the density profile is flat away from the pair of macrospheres (i.e. in what we can 
consider to be the bulk) and is equal to p. To keep p constant for each value of R, the number N is 
tuned to compensate the variation in the volume accessible to the small spheres. Periodic boundary 
conditions and the NVT ensemble are used. 

At the beginning of the simulation microparticles are placed randomly. At each MC step a mi- 
crosphere is selected at random and a random displacement is attempted. The displacement is 
accepted or rejected according to the Metropolis scheme, the displacement is tuned to reach an ac- 
ceptance ratio of 0.25 and cell lists are used in order to improve the efficiency of the algorithm. 30 
Configuration samples are taken after the mean square displacement of the microspheres equals 
O22 or after a number of moves sufficient to decorrelate the total energy per particle (in the SW and 
SS case). The value of the integrals that appear in Eq. (3) can be obtained using the relation 

2kv 2 p(S v ;R) sinO cos9 d9^(— £ cose,-), (18) 

•'° " r V<Vi<V+dv 
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where we sum all the cos ; whose distance v ; from the center of the macrosphere at position 
(R, 0, 0) is in the interval [v, v + dv] and the angle brackets denote an average over all the samples 
collected during a simulation. To obtain a better estimate of the integrals in Eq. (3) a third-order 
polynomial was fitted with the RHS of Eq. (18) corresponding to different v's and extrapolated the 
curves to v = a^ 2 and v = Ac^. The errors of the estimates are evaluated on the basis of the errors 
of the fit parameters. Data points in the fits were weighted according to the statistical uncertainties 
of the measured values of the RHS of Eq. (18). 



Plots of the density 

It's useful to obtain plots of the density obtained with the different methods in order to highlight 
any differences between their predictions. Such plots can be obtained from MC simulations by 
dividing the simulation box in cells, counting the number of particles contained in each of them, 
and averaging on multiple configurations. This data can be projected in 2D afterwards exploiting 
the symmetries of the system (for example p(r; R) =p(r, 0,0; R) = p(r,0 ; R) if defines a rotation 
around the direction of/?). 

Plots of the density associated to the integral equation and superposition method can be ob- 
tained using Eq. (13) by means of sampling it on a very fine mesh in real space. Down-sampling 
to the same mesh used with MC data produces results that can be compared to those obtained with 
the MC method. 

3.4 Numerical details 

We used the same geometries analyzed by Dickman et al. in 22, i.e. two size ratios E, = Oii/022 
= 5, 10, taking 022 to be of the form of Eq. (1) with 022 = 1 and £22 = 1 (that is, the smaller 
species is formed by unit hard spheres). The bulk packing fractions of the microspheres were also 
chosen in order to match those in (22): r\ = npo^/b = 0.1 16, 0.229, 0.341. Such choice of the 
parameters allowed us to test the validity of our data against (22). In addition we introduced a well 

12 



(shoulder) big-small interaction of the kind of Eq. (1), using A012 = 3.5, 5.5 (respectively in the 

case % = 5, 10) and whose depth (height) |e| = 1//3 = 1. 

As for the MC simulations, the size chosen for the box (H = 22, L = 16 for § = 5, H = 30, L = 24 

for § = 10), was big enough to keep bulk densities within the target values above with a precision 

of about 1% in all cases. Runs with different sizes and equilibration lenghts were performed in 

order to keep size and transient effects under 1%. 

Values of the reduced force f^ c {R) = PfMc(R)/(ftpOu) f° r values of R ranging from an to 

0\ i + 3 .0 022 at regular intervals of 0.2 022 were obtained. Each data point took about 10-20 hours 

of CPU time on an Intel Xeon 3.2 GHz processor to be determined. They are shown in Figure 6 

and Figure 7. Data taken from (22) (relative to the case e = 0) are also plotted for comparison. 

Table 1: Values of the parameters used to obtain the f^[ C (R) profiles. By a MC step here we mean 
a single particle displacement attempt. The averages have been evaluated by using the reported 
sampling frequency over the total number of MC steps shown in the last column. In all cases at 
least 10 8 equilibration MC steps were performed. 





N Sampling frequency (MC steps) Production (MC steps) 


£=5, tj =0.116 


« 1200 10 4 10 1U 


^ = 5, 77 = 0.229 


^2360 4-10 4 2-10 10 


£=5, 77 = 0.341 


^3520 5-10 5 4-10 10 


§ = 10, tj =0.116 


« 2700 10 4 10 10 


§ = 10, 77 = 0.229 


^5340 4-10 4 2-10 10 


§ = 10, 77 = 0.341 


^7950 5-10 5 4-10 10 



In the PY+S method values of gi2{ r ) were sampled on an equispaced (dr = 0.01) mesh of 4096 
or 8192 points in r-space respectively for § = 5, 10. Values for any r were obtained from the dis- 
crete sample through linear interpolation (linear extrapolation was used to obtain the values near 
discontinuous points). These in turn allowed us to obtain ppy+s( r ;^) f° r an y r through Eq. (13). 
The /py+s W cou ld finally be obtained performing a numerical integration of the RHS side of 
Eq. (3) for various values of R. The force profiles are shown in Figure 6, Figure 7. Each required 
a few minutes of CPU time on a desktop computer. The same procedure was carried out using the 
HNC+S method, and results are shown in Figure 6 for the case § =5, 77 = 0.1 16. 
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Forces profiles in the DL approximation are the dotted lines in the Figure 6 and Figure 7 and 
were obtained by straightforward substitution of the parameters above in Eq. (34). 

The phase diagram on the plane (t]i, rj) was obtained in the particular case of ^ = 5 for dif- 
ferent values of £12 and is shown in Figure 5(b). The associated effective potentials (shown in 
Figure 5(a)) were obtained by numerical integration of samples of the force in Eq. (34). The first 
order perturbation in in Eq. (10) was also found by numerical integration, while the derivative re- 
quired by in Eq. (11) was obtained via numerical differentiation. The discretization of real space 
needed to carry out such operations was performed on a grid sufficiently fine so that further refine- 
ments had no appreciable effect on the scale of the plots. 

In addition we obtained further information about the density of microspheres for the case 
t, = 5 and r\ = 0.116. 2D plots were obtained projecting the MC data in two dimensions exploit- 
ing the azimuthal symmetry of the problem. Equivalent diagrams were obtained via the PY+S 
method. In Figure 8, we show plots of [pMc(r,9;R) — Ppy+s(/ 5 0;^)]/Ppy+s(/ 5 0;^) at R = 5.2. 
Such diagrams allow to examine the differences between exact (though noisy) MC data and the ap- 
proximate PY+S predictions. Such plots are obtained setting pMc(r,0;R) = pMc(r,—0;R) when 
6 < 0, and noise at 6 rs is due to the bad statistics of the particles counts in this region (which is 
related to the small size of the bins = 2nr 2 sin drdO). 

4 Results 

The MC method allows to measure the exact effective force in the various cases. Our data are 
not always in good agreement with those found in (22) for the case £12 = and differences up 
to 20% are observed. Several simulations with different box sizes and production durations were 
performed and all our values obtained were consistent between each other within statistical error 
(which is way lower than 20%), confirming that our data are reliable. 
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In all cases examined here the DL approximation fails to describe MC data quantitatively. Still its 
results are in qualitative agreement with the exact profiles which show the features described in 
Section 3.1. Better results can be obtained via the PY+S approximation, that describes well the 
force profiles for £12 = and yields a reasonably good agreement also for £12 = ±1 at low densities 
and largest distances. 

Examination of the plots of the differences in densities with the MC and PY+S method elucidates 
the origin of the differences obtained in the force profiles. The satisfactory result obtained when 
£12 = is mirrored by a quite accurate match between Pmc and ppy+s- In this case PY+S only 
slightly underestimates the particle density in the zone close to both spheres (see Figure 8). The 
fact that the density match is very good everywhere but in this region suggests that this discrepancy 
is due to the superposition approximation. In the £12 = ±1 case the density differences are much 
more pronounced, again in the region close to both the macrospheres (indicating a breakdown of 
the superposition approximation) but also everywhere else in the vicinity of a single macrosphere 
(due to loss of accuracy of the PY closure). Similar results (shown here only in the case £ — 5, r\ = 
0. 1 16) were obtained via the HNC+S method, confirming that the closure plays a lesser important 
role than the superposition approximation in the disagreement with the MC force profiles. This 
last point is evident in the case r\ = 0. 1 16 of Figure 6, where the PY+S and HNC+S methods both 
fail to describe accurately the force profile at short ranges. 

The phase diagram (see Figure 5(b)) obtained from first order thermodynamic perturbation theory 
within the DL approximation shows that for small positive values of £12 one needs to move at 
higher densities of the smaller species as £12 increases in order to observe phase separation. This 
is due to the fact that in this regime an increase in £12 corresponds to an additional repulsive term 
in the effective potential (see Figure 5(a)). 

Such "stabilizing effect" due to mutual attraction doesn't hold at higher values of £12, where 
increments in £12 correspond to progressively stronger effective attraction and lowering of the 
coexistence line in the (rji, r\) plane. The critical density in the case under examination is higher 
in the low £12 regime than in the high £12 one. 



15 



5 Conclusions 

In this paper we have studied the effective forces between two big hard spheres dispersed into a 
fluid of smaller particles. In particular we studied the effect of interspecies interactions. To this 
aim we have extended the force derivation done by Attard to the case in which these take the shape 
of a square well or shoulder, that is in the case where both entropic and energetic effects play a 
role. As in the case of the original formula, it is sufficient to know the density profile of the small 
particles around the big ones to obtain the exact value of the effective force. The determination of 
the density however is a difficult task and one has to rely on some approximations. Here we have 
proposed an equivalent of the Asakura-Oosawa approximation, i.e. the assumption that the small 
particles behave like an ideal gas. This approximation leads to an analytical expression that cap- 
tures, at least qualitatively, the correct behavior as predicted by our MC simulations. An approach 
based on integral equations and a superposition approximation improves the results but stills fails 
in estimating the density in the neighborhood of the two large spheres regardless the choice of the 
closure used (at least in the cases examined here). We have shown how the addition of a shoulder 
or of a well in the interspecies interaction can change qualitatively the well known depletion force 
profile observed in hard sphere binary mixtures. The presence of a well, for example, causes the 
onset of a repulsion at short ranges and, in the case of deep wells, of an attraction at intermediate 
ranges. This happens because when a well is present the larger particles are surrounded by a layer 
of smaller ones. This layer makes a close contact between the large particles unfavorable, but at 
the same time if two bigger spheres share part of their surrounding layers the small particles sit- 
ting between them act as "glue", stabilizing a configuration of intermediate distance between the 
pair. Such description qualitatively agrees with what has been observed in molecular simulations 
of binary mixtures of eye-lens proteins 19 ' 20 of hard-core potentials with a Yukawa tail 31 and could 
be confirmed by experimental determinations of the force on colloids interacting with a square 
shoulder/potential (polymer grafted colloids might be a valid candidate member to this class) us- 
ing optical tweezers. 32 
The phase diagram as studied with simple thermodynamic perturbation theory within the DL ap- 
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proximation shows that the introduction of shallow well (low £12) has the effect of pushing to 
higher densities of the smaller component the phase separation. Increasing further £12, however, 
one reaches the point where deepening the well has the opposite effect. 

The effect of steric hindrance of the small particles in the purely hard sphere case has already 
been claimed to be used to stabilize solutions 33 and foods. 34 It's clear that taking into account the 
possibility of tuning the interspecies interactions could broaden even more the possible routes for 
stabilization. Summarizing, we confirm, in agreement with previous work, 19 ' 20 that the mutual 
attraction could be an extra parameter to play with when tuning the stability of a binary mixture. 
The present work provides a qualitative and quantitative analysis of the resulting changes. 
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A Derivation of the effective force 

It is shown in (21) that the force can be expressed as 

where 0n is the interaction potential between particles of type 1, F 2 is the free energy of the 
particles of type 2 that move in the potential generated by the fixed pair of type 1, Z 2 the partition 
function associated to it and /3 = 1/ksT. 

Following (21) we also have that 



dZa 
dR 



j^^-e- pV{r) )e^p(r)Z 2 dr, (20) 
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where V(r) is the potential energy of a single particle of type 2 due to the presence of the fixed 
couple of type 1. V(r) can be written as 



V(r) 



°o if r < On or \r — R\ < 0\2 

— £i2 if C12 < r < Xon xor On < \r — R\ < X<5\2 

-2ei2 if On. <r< Xo\2 and On < \r — R\ < ko\2 

otherwise 



(21) 



We can write the resulting Mayer function as: 

+ {e 2^2_ e Pen ) [^ an {r)^ Xon {r-R) + ^ Xan {r)^ an {r-R)] + 

-e 2 ^J^ 12 (r)J^ an (r-R)+ 

+ e^ [JT ai2 (r) + Jf ai2 (r-R)}, (22) 



where the geometric definition of the support of V (r) has been encoded using the characteristic 
functions J^b(f) defined as: 



Jtb(r) 



1 ifr<D 



ifr>D 



(23) 



Plugging Eq. (22) inside Eq. (20) and rearranging yields: 
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d -^ = -(l-e^)fl.^8(\r-R\-Xa l2) e^r) p{r)Zldr+ 

-(2e^-e 2 ^-l) fj^ 0l2 (r)^-^— ^ 8(\r-R\ -Xo n ) e^p{r)Z 2 dr+ 
J k \r — k\ 



_ (e 2/3 £l2 _ e /3 £l2) [jg> ( r) *.JL-± 8{\r-R\ -Xa l2 )e pv ^p(r)Z 2 dr+ 
J K \r — K\ 

_ {e 2Pen _ e^") J jr Xau (r)* • ^- 8(\r-R\ - o l2 )e^p{r)Z 2 dr+ 

+ e 2 ^ fjt^ir)^ ■ J— 1 8(\r-R\ - a l2 ) e" v ^p(r)Z 2 dr+ 
J K \r — K\ 

~ epEn [^■^ Z ^S(\r-R\-a i2 )eP v ^p(r)Z 2 dr, (24) 

J K \r — K\ 

where we have used the shorthand notation Xo\ 2 to refer to X\ 2 0\ 2 . 

The integrals in Eq. (24) can be evaluated in the three dimensional space performing the change 
of variables s = r — R, using spherical coordinates centered in R with 6 defined by s R = sRcos, 
and exploiting the azimuthal symmetry of the problem. 

The first one thus reads 

J R :- S -8(s-X<J l2 )e^ s+ Vp(s+R)Z 2 ds = 

= / / -■ -8(s- Xon) eP v{s+R) p(s+R)Z 2 s 2 sine dsd9d$ = 
Jo Jo Jo R s 

= 2nX 2 e 2 2 cosde pV{Sx °^p(S Xan )Z 2 sm6d6, (25) 

where in the last line Sx an is defined with the notation 

S v = v + R with|v|=v. (26) 

The sixth integral in Eq. (24) has the same form and its value is 

2ito 2 2 f cos6eP v{Sa i2)p(S ai2 )Z 2 sm6de. (27) 

Jo 
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As for the second 

R r-R 



j^ Xan {r)--*^5{\r-R\-Xon)e^p{r)Z 1 dr = 

2nX 2 o 2 2 f ' je Xan {S Xa JcoseeV v{s ^n)p(S Xai2 )Z 2 smede. (28) 



whereas the fourth 

R r-R 



J^Xo n (r)-.^- Wl 8(\r-R\-o i2 )eP v ^p(r)Z 2 dr = 
= 2ko1 2 j* ^ Xan (S Gn )coseeP y{ - s °n) p (S an )Z 2 smede . (29) 

while the third and the fifth ones vanish because wherever p is nonzero ffla xl is zero and viceversa. 
The ^"s can be eliminated introducing the auxiliary functions 



1 PACT1 2 (^C71 2 ) 



l pXo n {Sxo n ) 



P (s Xoi2 ) ifo<e< l e Xon 

otherwise 

p(s Xan ) ii l e Xan <e< 2 e Xai2 <K 

otherwise 



p(S ai2 ) ifO<0<X 2 
PanK^on) ~ \ > 

otherwise 



2 p an {S an )={ ~ , (30) 

otherwise 
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where 



(#-\ 



9 kai2 - arccos ^ -^— ^ 

&T,,, = arccos ' 



7 ^12 



0<t 12 = arccos 



-2AC712^ 

R 2 \ 



■ o n R) 



(3D 



are the angles at which the spheres of radii <j\2 an d h<3\2 centered in and R intersect. Defining 



p{Sxo n )= Pl<? n {Sxoi 2 )+ PXa n { S Xa l2 ) 

P\Son) = P(Ji 2 {S<T n ) + Po l2 {S<7 n ) 



(32) 



and substituting these inside the integrals and some algebra eventually yields Eq. (3). 



B Expression of the force in the dilute limit 



Use of Eq. (4) together with Eq. (21) allows to rewrite the densities that appear in Eq. (3) as 



P(SlcT 12 ) 



P(S 



on j 



p ifO<0< 1 A(Tl 

pe^ if l d Xan <d< 2 d x 

otherwise 

pe^ 1 ifO<0< 1 (7l2 

pe 2 ^ i{ i O an <O< 2 ai2 

otherwise 



0"12 



(33) 
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where the angles are the same defined in Eq. (31). Substitution of such density inside Eq. (3) yields 
a piecewise expression of the force valid for R e [<7n , +°°] : 



if an <R< Ion : 



/dl(*) 



f{e^U R2 + °\- X2 ^ ) + af 2] + 



+e 



2P £, 2 



R 



2 +1 / ^ + df 2 -A 2 of2 N2N 



' 12 

27? 



+ 



+(l-e /3£l2 : 



/? 



+e 



]8ei2 



+ A 2 a 1 2 2J + 

7? 2 + A 2 (y 2 2 _ (y 2 2 

27? 



+ 



/? 



if 2<7i2 <i? < (7i2+A(7i2 : 



/dl(R) 



-ff^(- f + 1; AV2 ) + ^i + 



+e 2 ^[-C7 2 2 + 



+(i-^ £ > 2 : 



i? 2 + a 2 9 -A 2 a 2 , 



'12 



12 



27? 



+ 



/? 



+e 



J8ei2 



+ A 2 C7 1 2 2 j + 

7? 2 + A 2 ct 2 2 -ct 2 2 
27? 



+ 



R 



if (7i2 + A(7i2 <i? < 2A<7i2 : 



fyL(R) = -JUl-ef>^ 



R 



+AV? + 



12 



+^(-^+(£f)]} 
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if2A<7i2</g| : 

Idl(R)=0. (34) 



C Linear dependence of volume terms with p 



1 



Following Dijkstra et al. 23 we calculate the volume terms (0- and 1-body) of Q. in the effective 
potential of Eq. (6). 

The first term can be interpreted as the grand potential of a pure system of small particles at fugacity 
Z2 enclosed in a volume V 

^=f// r =f' (35) 

while the second is 

= lf(y^((^ £l2 -i)(^ 3 -i)-i)) (37) 

= piV— x const, (38) 

where we have used the definition ft = e^^ l2 ^ Ri ~ r > — 1. From Eq. (35) and Eq. (38) we see that 
(Q.q + Q.\)/V is linear with respect to p\ and thus does not alter the results of the phase diagram 
construction that leads to Figure 5(a). 12 ' 23 
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Figure 1: Scheme of the geometrical parameters that appear in Eq. (2) and Eq. (3). The solid 
white area indicates the region inaccessible to the centers of the microparticles due to the hard part 
of the potential 012(f)- The rippled area represents the structuration of the density p(r) due to the 
presence of the macrospheres. 
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Figure 2: Qualitative force profiles for various values of £12 in the dilute limit. This qualitative 
behavior is observed in a quite broad region in the space of the geometrical parameters. The 
configurations corresponding to R = 0\\,2<Ji2,<5\2 + Xo\2,2Xo\2 are shown. The smaller disks 
represent the hard spheres, the white coronas the volume forbidden to the centers of the micro- 
spheres, and the bigger coronas are the square well/shoulder. The size of a microsphere is that of 
the red disks. The geometrical parameters used in the drawing are an = 5, O12 = 3, Xon = 3.75. 
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(a) 



(b) 




(c) 



(d) 



Figure 3: Figure 3(a), Figure 3(b): Plot of the density of small particles around a big sphere (the 
blue disk) at large R (isolated particle case) and when another sphere (not shown) is present at 
small R in the case £12 = and in the dilute limit. The anisotropy of the density determines the 
onset of an attractive force. 

Figure 3(c), Figure 3(d): Density of small particles around a macrosphere at large R and small R 
in the case £12 < 0. The anisotropy of the density determines the onset of an attractive force as in 
Figure 3(b), this time due to the contribution of the density both at the surfaces of radius Ou (hard 
core) and Xou (outer rim of the shoulder). 
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(a) 



(b) 




(c) 



(d) 



Figure 4: Density of small particles around a macrosphere for decreasing R in the case £12 > 0. 
The effective interaction is the result of the interplay between the "push" due to the density at the 
surface of radius On (much stronger due to the higher density inside the well) and the "pull" at 
A(7i2 (outer rim of the well). 



30 



f = 5 




(a) 



0.45 



0.40 



0.35 



0.30- 



0.25- 



0.20- 



0.15 



\jF 


+ + + + 




'" T T T" 


\y " 


-' 


. . £i 2 = 0.0 


- 


t t e 12 = 0.05 - 




£12 = 0.1 




□ □ £ i2 = 1.05 - 




o o £i2 = 1.1 


- 


* * e 12 = 1-15 - 


, 



0.1 



0.2 0.3 



0.4 



(b) 

Figure 5: Figure 5(a): Effective potentials obtained by numerical integration of the force in Eq. (34) 
using the parameters reported in the text relative to the case ^ = 5 for various values of £12. Fig- 
ure 5(b): Coexistence curves obtained by first order perturbation theory using the potentials in 
Figure 5(a). 
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Figure 6: Force profiles obtained with the various methods in the case E, = 5, 7} = 
0.116, 0.229, 0.341 for different values of £12- In the case 77 =0.116 the force profiles obtained 
with the HNC closure and the superposition approximation are also presented. The legend that 
holds for all the other force profiles in the article is that shown for the case r\ = 0.229. 
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Figure 7: Same as Figure 6, but using t, = 10. MC data are obtained using a simulation box whose 
H = 30 and L = 24. 
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Figure 8: Plot of the difference of the densities of smaller particles around a pair of large particles 
obtained via MC and PY+S for different values of £12. 
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